An Improved Numerical Integration Method for Flight Simulation

In this paper a modified form of Euler integration is described which, when applied to the six-degree of freedom flight equations, retains and enhances many of the advantages of AB-2 integration and at the same time eliminates the disadvantages. The scheme is based on the Euler integration formula, but with the state-variable derivative represented at the midpoint of each integration step. In this case the conventional first-order Euler method actually becomes second order, with a very small accompanying error coefficient. To apply this method to the six-degree-of-freedom flight equations it is necessary to define velocity states at half-integer frame times and position states at integer frame times. It is shown through dynamic error analysis that the modified Euler method has an error coefficient which is one-tenth that associated with AB-2. The method also exhibits minimal output delay in response to transient inputs. The modified Euler method may also be useful in the integration of state and costate equations in real-time mechanization of Kalrnan filters for navigation and control systems. The ever increasing complexity of the math models used as a basis for real time flight simulation has continued to apply pressure on digital processor speed requirements for such simulations. More effective numerical integration algorithms can help relieve some of this pressure. The most popular method currently in use for flight simulation is the AdamsBashforth second-order predictor method, usually referred to as AB-2. Its advantages include second-order accuracy with respec1 to integration step size, only one required pass through the state equations per integration step, and compatibility with real-time inputs. Disadvantages include stability problems associated with extraneous roots and response delays of one or two frames following transient inputs. In the next section we consider AB-2 along with two other second-order integration methods suitable for real-time simulation of dynamic systems. The first is a two-pass realtime predictor-corrector algorithm and the second is a singlepass version of the same. We then introduce the modifiedEuler method and describe its application to the flight equations. The accuracy of the various methods is compared by means of time-history plots of the dynamic error in simulating aircraft response to a control-surface input. 2. Some Real-time Second-order Inteeration Aleorithm~ Consider first the AB-2 predictor integration method applied to the state equation given by Professor of Aerospace Engineering Associate Fellow, AIAA Copyright 0 American Institute of Aeronauticd wd Astronautics, Inc., 1989. AU rights reserved. Here X is the state vector, U is the input vector, and h is the integration step size. The standard AB-2 predictor algorithm is the following: where X , = X(nh) and The AB-2 formula in Eq. (2) is derived from the area under a linear extrapolation from F, to Fn+l based on Fn and Fn.l. From Z transform theory it can be shown that the numerical integration formula in general has a transfer function for sinusoidal inputs which takes the form [I] where for AB-2 integration, k = 2 and el = 5/12. The term eI(iwh)k represents the error in the integrator transfer function compared with the ideal continuous integrator transfer function, 1 Based on the integrator model of Eq. (4) the transfer function gain and phase error in simulating any order of dynamic system, when quasi-linearized, can easily be obtained [I]. It can also be shown that the fractional error in any characteristic root as a result of integration truncation errors is given by where 1 is the continuous system root and a* is the equivalent root for the digital simulation. Thus it is apparent that the order k and error coefficient eI for any given integration algorithm can be used to pred~ct he dynamic errors that will be i n d u c e d into a simulation because of the finite step size h when using that algorithm. It should be noted that this methodology of error analysis is not applicable to multiple-pass integration methods, such as the Runge-Kutta algorithms, where different orders of integration algorithms are used in the various derivative evaluation passes that constitute a single integration step. It is applicable to all of the real-time methods considered in this paper. As indicated in Eq. (1). the state variable derivative F will in general depend on the state X. Since the AB-2 algorithm in Eq. (2) involves the past derivative Fn.1, the next state Xn+1 will depend on the past state Xn-1 as well as the current state X n . For this reason the AB-2 method introduces one extraneous state per integration. The characteristic roots corresponding to these extraneous states damp rapidly for small integration step sizes, in which case they do not contribute significant errors to the simulation. However, when the step size becomes large, the extraneous roots can cause instability. In particular, for a negative real root A, instability occurs when Ah < -1. This means that when AB-2 integration is used, the step size must be kept less than the shortest time constant in the system being simulated. Stability charts in the complex ah plane show the allowable step sizes when the roots of the continuous system are complex [2]. Consider next the Adams-Moulton two-pass predictorcorrector algorithm. Here the AB-2 method is used on the first pass to compute an estimate, &+I, for the n+l state. From this state and the input Un+l the estimated derivative is in turn calculated. The corrector pass then computes Xn+1 using Fn and with the following formula: If the estimate in Eq. (6) is repaced by the true derivative F,+I, the formula represents implicit trapezoidal integration. It turns out that the explicit AM-2 method has the same asymptotic error coefficient, e, = 1/12, as the implicit trapezoidal integration which it approximates. The order k = 2. It should be noted that the AM-2 algorithm requires two passes through the state equations per integration step. At the beginning of the second pass, the input Un+1 is used in the calculation of the derivative estimate FA+1. But Un+1 will not be available in real time until the completion of the second pass. Hence AM-2 is not compatible with real-time inputs. However, a second-order predictor-corrector algorithm suitable for real-time inputs can be constructed using the concept behind modified Euler integration. In modified Euler integration the state-variable derivative is represented at the midpoint of the integration step. Thus the integration formula becomes Xn+l = Xn + h h+ln (7) If the derivative F is a function of the state X, which is normally the case, then it is necessary to compute an estimate for the state Xn+1n in order to evaluate Fn+lR. In real-time Runge-Kutta 2 this is accomplished using Euler integration with a step size of h/2. The resulting real-time RK-2 integrator error coefficient el = 116. In the second-order real-time predictor-corrector method, denoted here as RTAM-2, the estimate forXn+1n is computed using a second-order predictor algorithm. This leads directly to the following difference equations [31: x&,, = x, + h(QFn-+Fn, ) (8) The predictor formula in Eq. (8) for Xn+1n is derived from the area under a linear extrapolation based on Fn and Fn-1. The RTAM-2 given by Eqs. (8), (9) and (10) requires the input Un at the start of the first pass and U,,+ln at the start of the second pass, both compatible with real time. Here the RTAM-2 integrator error coefficient el = 1/24, compared with -1112 for standard AM-2. For a negative real root A, instability occurs when Ah < 2. In the complex ah plane the overall stability boundary is slightly larger than the stability boundary for standard AM-2. Thus the modified AM-2 integration denoted here as RTAM-2 can be used as a general, real-time integration method for simulating nonlinear systems. For small step sizes h it is twice as accurate as either implicit trapezoidal integration or standard AM-2 integration, neither of which are compatible with real-time inputs. It is ten times more accurate than AB-2 integration. It should be realized, however, that the RTAM-2 considered here is a two pass per step method. Thus it will normally take about twice as long for computer execution as AB-2. When this speed differential is taken into account, the modified AM-2 still exhibits 2.5 times the dynamic accuracy of AB-2 based on the approximate asymptotic formulas for small step size. In a real-time simula-tion the intermediate state X;1+112 can be used as a real-time output which has full secondorder predictor accuracy for the step size h/2. Use of both X;1+in and Xn+1, then, provides outputs at the sample rate of a single-pass method, even though a two-pass method is being utilized. Note also that the method uses two input samples per frame, Un and Un+1n. The final integration method considered in this section is a single-pass version of the RTAM-2 algorithm described above. The method computes the state-variable derivative F only at integer frame times rather than both half-integer and integer frame times, as in Eqs. (8), (9) and (10). Values of the state X at half-integer frame times are computed using the modifiedEuler algorithm, while values of X at integer frame times are computed using a second-order predictor. The difference equations are the following [3]: